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Abstract 

The non-linear response of infinite periodic solids to homogenous electric fields and collective 
atomic displacements is discussed in the framework of density functional perturbation theory. The 
approach is based on the 2n+l theorem applied to an electric- field-dependent energy functional. We 
report the expressions for the calculation of the non-linear optical susceptibilities, Raman scattering 
efficiencies and electrooptic coefficients. Different formulations of third-order energy derivatives are 
examined and their convergence with respect to the k-point sampling is discussed. We apply our 
method to a few simple cases and compare our results to those obtained with distinct techniques. 
Finally, we discuss the effect of a scissors correction on the EO coefficients and non-linear optical 
susceptibilities. 

PACS numbers: 78.20.-e,71.15.Mb,78.20.Jq 
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I. INTRODUCTION 



Nowadays, density functional theory^'^ (DFT) is considered as a standard method in 
condensed matter physics, to study electronic, structural and macroscopic properties of 
solids from first principles. Combined with adiabatic perturbation theory, it allows a priori 
the computation of derivatives of the energy and related thermodynamic potentials up to 
any order. At the second order, this approach has been applied to compute linear response 
functions such as phonon frequencies or Born effective charges with an accuracy of a few 
percents. The third-order derivatives are related to non-linear properties such as phonon 
lifetimes, Raman tensors or non-linear optical susceptibilities. 

The linear-response formalism has been implemented in various first-principles codes and 
is routinely applied to an increasing number of systems [see for example Ref. [3] and ref- 
erences therein]. By contrast, the non-linear response formalism has been mostly restricted 
to quantum chemistry problems. Although the hyperpolarizabilities of a huge number of 
molecules have been computed, taking into account both electronic and vibrational (ionic) 
contributions^'^, applications in condensed matter physics have focused on rather simple 
cases^"-^^. 

Here, we present a methodology for the computation of nonlinear response functions and 
related physical quantities of periodic solids from density functional perturbation theory 
(DFPT). We focus on perturbations characterized by a zero wavevector and involving either 
three electric fields, or two electric fields and one atomic displacement. Following Nunes and 
Gonze^^, our approach makes use of the 2n+l theorem applied to an electric- field-dependent 
energy functionaP^. We report the local density approximation (LDA) expressions, as im- 
plemented within the ABINIT package^^. 

Our paper is organized as follows. In Sec. II, we describe the theoretical background 
related to the 2n+ 1 theorem and the electric field perturbation. In Sec. Ill, we describe the 
computation of the non-linear optical susceptibilities, the non-resonant Raman scattering 
efficiencies of both transverse (TO) and longitudinal (LO) zone-center optical phonons and 
the hnear electrooptic (EO) tensor. In Sec. IV, we illustrate the validity of the formal- 
ism by applying our methodology to some semiconductors and ferroelectric oxides and we 
briefly discuss the effect of a scissors correction on the EO coefficients and non-linear optical 
susceptibilities. 
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Some of the tensors we consider in this work depend on static electric fields : they include 
contributions of both the electrons and the ions. Other quantities imply only the response of 
the valence electrons: they are defined for frequencies of the electric field high enough to get 
rid of the ionic contributions but sufficiently low to avoid electronic excitations. For clarity, 
we adopt the following convention. Static fields will be labeled by greek indices {a, P, ...) 
while we will refer to optical fields with roman symbols (i, j, ...). To simplify the notations, 
we will also drop labels such as " oo" for quantities that do not involve the response of the 
ions. Using this convention, we can write £y and s^/?) respectively, for the optical and static 
dielectric tensor, respectively, and Vij^ for the linear EO (Pockels) tensor that involves two 
optical and one static electric field. 



II. FORMALISM 



A. Mixed third-order energy derivatives 



In this section, we present the general framework of the computation of third order energy 
derivatives based on the 2n+l theorem^'''"^^. Using the notations of Refs. [20,21], we consider 
three Hermitian perturbations labeled Ai, A2 and A3. The mixed third-order derivatives 
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T is the kinetic energy and Ehxc (resp. vhxc) is the sum of the Hartree and exchange- 
correlation energy (resp. potential). The first-order potential v^^^ can be computed as a 
second-order functional derivative of Ehx(^^ 



d\2 5n{r) 



(4) 



A=0 



Within the parallel gauge, the first-order Lagrange multipliers are given by 

As a consequence of the "2n + 1" theorem, the evaluation of Eq. (3) requires no higher 
order derivative of the wavef unctions than the first one. These first-order wavef unctions 
are nowadays available in several first-principles codes. They can be computed from linear 
response by minimizing a stationary expression of the second-order energy^° or equivalenty 
by solving the corresponding Sternheimer equation^^. It follows that the computation of 
third-order energies does not require additional quantities than the calculation of second- 
order energy derivatives. 

Eq. (3) is the general expression of the third-order energy derivatives. In case at least 
one of the perturbations does not affect the explicit form of the kinetic energy or the Hartree 
and exchange-correlation energy, it can be simplified : some of the terms may be zero. This 
is the case for the electric field perturbations treated in this work as well as for phonon 
type perturbations. Further simplifications can be made in case pseudopotentials without 
non-linear exchange-correlation core-correction are used. 



B. The electric field perturbation 

As mentioned in the introduction, special care is required in case one of the perturbations 
\j is a macroscopic electric field £. In fact, as discussed in the litterature, for infinite periodic 
solids, usually treated with Born- von Karman boundary conditions, the scalar potential £ ■ r 
breaks the periodicity of the crystal lattice. Moreover, it is unbound from below : it is 
always possible to lower the energy by transferring electrons from the valence states to the 
conduction states in a distant region (Zener breakdown). However, for sufficiently small 
fields, the tunneling current through the band gap can be neglected and the system is well 
described by a set of electric field dependent Wannier functions W„(r). As shown by Nunes 
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and Vanderbilt^^, these Wannier functions minimize the energy functional 

E [Wn\ S\ = Eo [Wn] -nS V (6) 

where Eq is the Kohn-Sham energy under zero field and the macroscopic polarization 
that can be computed from the Wannier function centers. It is important to note that these 
Wannier functions do not correspond to the true ground-state of the system but rather to a 
long lived metastable state. 

In practical applications, it is not mandatory to evaluate the functional Eq. (6) in a 
Wannier basis. It can cquivalcntly be expressed using Bloch functions Unk related to Wn by 
a unitary transform. In this case, the polarization can be computed as a Berry phase of the 
occupied bands^^ 

occ « 

\3 / C?k(li„k|Vfe|link) (7) 

„ Jbz 



^ 2ie 



{2ny 

where BZ is the Brillouin zone, e the absolute value of the electronic charge and the factor 
of 2 accounts for spin degeneracy. The Bloch functions are chosen to satisfy the periodic 
gauge condition 

e'^"^Mnk+G = Unk. (8) 

In order to use Eq. (7) in practical calculations, the integration over the BZ, as well as the 
differentiation with respect to k, have to be performed on a discrete mesh of k-points. In 
case of the ground-state polarization, the standard approach is to build strings of k-points 
parallel to a vector of the reciprocal space Gy. The polarization can then be computed 
as a string- averaged Berry phase. Unfortunately, the adaptation of this method to the 
computation of the energy derivatives is plagued with several difficulties, like the following. 
The general form of the non-linear optical susceptibility tensor of a compound is imposed 
by its symmetry. For example, in zinc-blende semiconductors, this tensor, expressed in 
cartesian coordinates reduces to Xifi — X^^^kyd where e is the Levi-Civita tensor. It follows 
that the reduced coordinates formulation of x^^] satisfy the relation 



(2) 



A2) 



where at least one of the three indices i, j and / are different from the two others. When 

('2) 

we tried to use strings of k-points to compute Xijii Eq. (9) was not satisfied. However, 
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we were able to avoid these problems, by using the finite difference formula of Marzari and 
Vanderbilt^^ on a regular grid of special k-points (instead of strings) 

V/(k) = Yl [/(k + b) - /(k)] (10) 

b 

where b is a vector connecting a k-point to one of its nearest neighbours and w^, is a weight 
factor. The sum in Eq. (10) includes as many shells of first neighbours as necessary to 
satisfy the condition 

$^«^bM/3 = -^ (11) 

where ba are the reduced coordinates of b and Qa/s is the metric tensor associated to the real 
space crystal lattice. 

In case of the ground-state polarization, we cannot apply the discretization Eq. (10) 
directly to Eq. (7). As shown by Marzari and Vanderbilt, the discretization of Eq. (7) does 
not transform correctly under the gauge transformation 

unk{r) ^ e-^'^-'^Unkir). (12) 

Since Eq. (12) is equivalent to a shift of the origin by one lattice vector R, V must change 
accordingly by one polarization quantum. In order to obtain a discrete expression that 
matches this requirement, we must combine Eq. (10) with the King-Smith and Vanderbilt 
formula^^'^^ 

^ = ]^ E E ^bb$> In det [5(k, k + b)] (13) 
k b 

where S is the overlap matrix between Bloch functions at k and k -|- b: 

S'n,m(k, k + b) = (M„k|Mm,k+b)- (14) 

As discussed by Nunes and Gonze^^, when we apply the perturbation expansion of the 
preceeding section to the energy functional Eq. (6), we can adopt two equivalent approaches. 
The first possibility is the use of Eq. (7) for the polarization and a discretization after hav- 
ing performed the perturbation expansion. The second possibility is to apply the 2n + 1 
theorem directly to Eq. (13) in which case no additional discretization is needed. Using the 
notations of Nunes and Gonze, we will refer to the first case as the DAPE (discretization 
after perturbation expansion) formulation and to the second one as the PEAD (perturba- 
tion expansion after discretization) formulation of the third-order energy. In the following 
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sections, we will discuss both expressions. In addition, in Sec. IV B, we will compare their 
convergence with respect to the k-point sampling on a realistic example. The perturbation 
expansion of the first term {Eq) of Eq. (6) can easily be performed as it is described in the 
Sec. II A. At the opposite, the expansion of the second term {—Q£ ■ T*) is more tricky since 
it explicitely depends on the polarization. In the two sections that follow, we will focus on 
the —0,£ ■ V term of Eq. (6). It will be referred to as Epoi. 

C. DAPE expression 

According to the formalism of the preceeding section, the £ ■ V term acts as an additional 
external potential that has to be added to the ionic one. The £ ■ V perturbation is linear 
in the electric field and does not depend explicitely on other variables such as the atomic 
positions. It just enters the terms of Eq. (3) that involve the first derivative of Vf,xt with 
respect to £. In other words, the only terms in Eq. (2) that involve the expansion of P are 
of the form E^^^^^^ where Ai and A3 represent an arbitrary perturbation such as an electric 
field or an atomic displacement. 

The DAPE expression of the third-order derivative of Epoi writes as follows 

where u^^^, are the projection of the first-order wavefunctions on the conduction bands. 
The complete expression of various third-order energy derivatives, taking into account the 
expansion of both Eq and Epoi, are reported in Sec. III. Eq. (15) was derived first by 
Dal Corso and Mauri^^ in a slightly different context: they performed the perturbation 
expansion of the energy functional Eq. (6) using a Wannier basis. The resulting expression 
of the third-order energy was expressed in terms of Bloch functions by applying a unitary 
transform to the Wannier orbitals. 

Using the finite difference expression of Marzari and Vanderbilt Eq. (10), Eq. (15) 
becomes 

fc 6 n,m 

- {<kAK,m} (16) 
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where Gi is a basis vector of the reciprocal lattice. 



D. PEAD expression 



Applying directly the 2n + 1 theorem to Eq. (13) we obtain the alternative PEAD 
formulation of the third-order energy: 
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where Q is the inverse of the overlap matrix 5" and 5"^^ its first-order perturbation expansion 



S^lUk, k + b)^ (u'r^Ju^X,) + (u'X^^,). 



(18) 



III. NON-LINEAR OPTICAL PROPERTIES 

In the previous section we have discussed the general expressions of third energy deriva- 
tives. We now particularize them to the computation of selected non-linear properties. 

A. Non-linear optical susceptibilities 



In an insulator the polarization can be expressed as a Taylor expansion of the macroscopic 
electric field 

^=^+Exs^^+i:4fe^^+--- (19) 

j=l j,k=l 

where P/ is the zero-field (spontaneous) polarization, Xif the linear dielectric susceptibility 
(second rank tensor) and Xi/i the second-order non-linear optical susceptibility (third rank 
tensor). In the litter ature on non-linear optics, one often finds another definition of the non- 
linear optical susceptibility : instead of x^^]-> it is more convenient to rely on the d-tensor 
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defined as 

diji^^X^l (20) 
In general, the polarization depends on valence electrons as well as ions. In the present 
section, we deal only with the electronic contribution : we will consider the ionic cores 
as clamped to their equilibrum positions. This constraint will be relaxed in the following 
sections where we allow for ionic displacements. 

Experimentally, the electronic contribution to the linear and non-linear susceptibilities 
corresponds to measurements for electric fields at frequencies high enough to get rid of 
the ionic relaxation but low enough to avoid electronic excitations. In case of the second- 
order susceptibilities, this constraint implies that both the frequency of £, and its second 
harmonic, are lower than the fundamental absorption gap. 

The general expression of the electronic non-linear optical susceptibility depends on the 
frequencies of the optical electric fields [see for example Ref. [27]]. In the present context 
of the 2n+ 1 theorem apphed within the LDA to (static) DFT, we neglect the dispersion 
of Xiji*- As a consequence, xifi satisfies Kleinman's^^ symmetry condition which means that 
it is symmetric under a permutation of i, j and /. In order to be able to investigate its 
frequency dependence, one would need either to apply the formalism of time-dependent 
DFT^ or to use expressions that involve sums over excited states^^"^^. Following the work 
of Dal Corso and co-workers^'^^ we can relate the non-hnear optical susceptibihties to a 
third-order derivative of the energy with respect to an electric field 

- -^E^-'^'^ (21) 
where E^^^^^'- is defined as the sum over the permutations of the three perturbations E^*^^^' 
(2). Using the PEAD formulation of Sec. II B we can express these terms as follows 
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B. Raman susceptibilities of zone-center optical phonons 



We now consider the computation of Raman scattering efficiencies of zone-center optical 
phonons. In the hmit qr — > 0, the matrix of interatomic force constants C can be expressed 
as the sum of an analytical part and a non- analytical term^^ 



The analytical part corresponds to the second-order derivative of the energy with respect 
to an atomic displacement at q = under the condition of vanishing macroscopic electric 
field. The second term is due to the long-range electrostatic interactions in polar crystals. 
It is at the origin of the so-called LO-TO sphtting and can be computed from the knowledge 
of the Born effective charges Z*^^ and the electronic dielectric tensor^^ Sij. The phonon 
frequencies Um and eigendisplacements Um{tiOi) are solution of the following generalized 
eigenvalue problem 



k',/3 

where is the mass of atom k. As a convention, we choose the eigendisplacements to be 
normalized as 



In what follows we consider non-resonant Raman scattering where an incoming photon of 
frequency ujq and polarization eo is scattered to an outgoing photon of frequency {ujq — cUm) 
and polarization 65 by creating a phonon of frequency cUm (Stokes process). The scattering 
efficiency^^'^^ (cgs units) corresponds to 



C^aMQ - 0) = Ct%{q = 0) + C,%(q ^ 0). 



(23) 




(24) 




(25) 




(26) 



where c is the speed of light in vacuum and rim the Boson factor 



1 



(27) 



n, 



The Raman susceptibility a™ is defined as 




(28) 
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where xlj is the electronic hnear dielectric susceptibihty tensor. Q is the angle of collection 
in which the outgoing photon is scattered. Due to Snell's law, Q is modified at the interface 
between the sample and the surrounding medium. Experimentally, the scattering efficiencies 
are measured with respect to the solid angle of the medium while Eq. (26) refers to the 
solid angle inside the sample. In order to relate theory and experiment, one has to take 
into account the different refractive indices of the sample and medium. For example, in the 
isotropic case, Eq. (26) has to be multiplied^^ by {n' /nf where n and n' are respectively 
the refractive indices of the sample and the medium. 

For pure transverse optical phonons, -q;^ can be computed as a mixed third-order deriva- 
tive of the energy with respect to an electric field, twice, and to an atomic displacement 
under the condition of zero electric field 



In case of longitudinal optical phonons with wavevector qr ^ in a polar crystal, Eq. (28) 
must take into account the effect of the macroscopic electric field generated by the lattice 
polar vibration. This field enters the computation of the Raman susceptibilities at two 
levels. On one hand, it gives rise to the non- analytical part of the matrix of interatomic 
force constants Eq. (23) that modifies the frequencies and eigenvectors with respect to pure 
transverse phonons. On the other hand, the electric field induces an additional change in 
the dielectric susceptibility tensor related to the non-linear optical coefficients xijk- ^'^'^ 
longitudinal optical phonons, Eq. (29) has to be modified as follows^^: 



The mixed third-order derivatives (29) can be computed from various techniques including 
finite differences of the dielectric tensor^'''"^^ or the second derivative of the electronic density 
matrix^°'^^. Here, we follow an approach similar to Deinzer and Strauch^° based on the 2?t,-|-1 
theorem. The third-order energy can be computed as the sum over the 6 permutations Eq. 
(2) of Tka, £i and £j. According to the discussion of Sec. II B, we have to distinguish between 
the terms that involve the discretization of the polarization such as E'^"'^^*^^ or E^^^^'^'^^ and 
those that can be computed from a straightforward application of the 2n + 1 theorem such 
as E^i'^'^^'^i. The former ones show an electric field as second perturbation. They can be 
computed from an expression analogous to Eq. (22) 




(29) 




(30) 
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(31) 



5n(r)(5n(r')5n(r") 

We obtain a similar expression for E^^^^'^'^^. The remaining terms do not require any differ- 
entiation with respect to k. They can be computed from the first-order change of the ionic 
(pseudo) potential with respect to an atomic displacement v1^^ 
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d S^E 



Hxc 



dr^x 6n{r)Sn{r') 



n^Hr)n^^(r') 



drdr'dr" 



\n 



—n^'^^{r)n^'{r')n^^{r"). 



(32) 



6n{r)6n{r')6n{r") 

In pseudopotential calculations, the computation of the first-order ionic potential vl^^ re- 
quires the derivative of local and non-local (usually separable) operators. These operations 
can be performed easily without any additional workload as described in Ref . [20] . 

In spite of the similarities between Eqs. (31) and (32) and the expression proposed 
by Deinzer and Strauch wc can quote few differences. First, our expression of the third- 
order energy makes use of the PEAD fomulation for the expansion of the polarization. 
Moreover, Eq. (32) is more general since it allows the use of pseudopotentials with non- 
linear core correction through the derivative of the second-order exchange-correlation energy 
with respect to t^x (third term). 
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C. Sum rule 



As in the cases of the Born effective charges and of the dynamical matrix^^, the coefficients 
dXif /dr^ce must vanish when they are summed over all atoms in the unit cell. 

E^ = (33) 

Physically, this sum rule guarantees the fact that the macroscopic dielectric susceptibility 
remains invariant under a rigid translation of the crystal. In practical calculations, it is not 
always satisfied although the violation is generally less severe than in case of C or Z*. Even 
in calculations that present a low degree of convergence, the deviations from this law can be 
quite weak. They can be corrected using similar techniques as in case of the Born effective 
charges^^. For example, we can define the mean excess of dx\f /dr^a per atom 



9x!ff _ 1 V- dx^ 



dTa Nat ^ dr, 
and redistribute it equally between the atoms 



= J-yp- (34) 



9y(^) dx^'^ dx^'^ 

D. Electrooptic tensor 

The optical properties of a compound usually depend on external parameters such as the 
temperature, electric fields or mechanical constraints (stress, strain). In the present section 
we consider the variations of the refractive index induced by a static or low- frequency electric 
field S^. At the linear order, these variations are described by the linear EO coefficients 
(Pockels effect) 

3 
7=1 

where {£~^)ij is the inverse of the electronic dielectric tensor and rjj-y the EO tensor. 

Within the Born and Oppenheimer approximation, the EO tensor can be expressed as 
the sum of three contributions: a bare electronic part rfj^, an ionic contribution rjj" and a 
piezoelectric contribution rfj^^°. 
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The electronic part is due to an interaction of £-y with the valence electrons when con- 
sidering the ions artificially as clamped at their equilibrium positions. It can be computed 
from the non-linear optical coefficients. As can be be seen from Eq. (19), xfji defines the 
second-order change of the induced polarization with respect to Taking the derivative 



of Eq. (19), we also see that x\jl defines the first-order change of the linear dielectric sus- 
ceptibility, which is equal to ^Ae^j. Since the EO tensor depends on A(£~^)jj rather than 
Asjj , we have to transform Asj^ by the inverse of the zero field electronic dielectric tensor^^ 

3 

m,n=l 

Using Eq. (37) we obtain the following expression for the electronic EO tensor 

4. = E i^~"hx'uk^~")i'3 , ■ (38) 
l,l'=l ' 

Eq. (38) takes a simpler form when expressed in the principal axes of the crystal under 

investigation^^ 

(39) 



„e; ^ -8^, ,(2) 

^ 3 



where the rii coefficients are the principal refractive indices. 

The origin of the ionic contribution to the EO tensor is the relaxation of the atomic 
positions due to the applied electric field £^ and the variations of Sij induced by these 
displacements. It can be computed from the Born effective charges ^* and the |^ 
coefficients introduced in Sec. IIIB. As shown in appendix A [see also Refs. 36,45], the 
ionic EO tensor can be computed as a sum over the transverse optic phonon modes at q = 



m 



where ol^ is the Raman susceptibihty of mode m [Eq. (28)] and Pm,p the mode polarity 

Pm,7 = E ^K,lP'^rn{l^(^)- (41) 

that is directly linked to the modes oscillator strength 

Sm,a(} = Pm,a " Prra,/3- (42) 

For simplicity, we have expressed Eq. (40) in the principal axes while a more general 
expression can be derived from Eq. (37). 
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Finally, the piezoelectric contribution is due to a relaxation of the unit cell shape due to 
the converse piezoelectric effect As it is discussed in appendix A, it can be computed 
from the elasto-optic coefficients Pijfj,i, and the piezoelectric strain coefficients (i^^,^ 

3 

^ijl ~ Pijfiud"ffiu- (43) 

/i,i/=i 

In the discussion of the EO effect, we have to specify whether we are dealing with strain- 
free (clamped) or strcss-frce (undamped) mechanical boundary conditions. The clamped 
EO tensor r^j^ takes into account the electronic and ionic contributions but neglects any 
modification of the unit cell shape due to the converse piezoelectric effect^^'^^ 

r'f. =r^' +r^°". (44) 

' ij7 iJ7 tJI \ / 

Experimentally, it can be measured for frequencies of high enough to eliminate the re- 
laxations of the crystal lattice but low enough to avoid excitations of optical phonon modes 
(usually above ~ 100 MHz). To compute the undamped EO tensor r^j^, we have to add the 
piezoelectric contribution to r^-^ 

Experimentally, r'^j^ can be measured for frequencies of £^ below the (geometry dependent) 
mechanical body resonances of the sample^^ (usually below ~ 1 MHz). 

IV. RESULTS 

A. Technical details 

Our calculations have been performed within the local density approximation (LDA) 
to the density functional theory^'^ (DFT). We used the ABINIT^^ package, a planewave, 
pseudopotcntial DFT code^^ in which wc have implemented the formalism presented above. 
For reasons that will become obvious below, we chose the PEAD formulation Eq. (17) 
to perform the differentiation with respect to k. For the exchange-correlation energy E^c 
we relied on the parametrization of Perdew and Wang^^ as well as the parametrization 
of Goedecker, Teter and Hutter^*^. These expressions have the advantage to avoid any 
discontinuities in the functional derivative of E^c- 
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In case of the semiconductors Si, AlAs and AlP, we used at 16 x 16 x 16 grid of special 
k-points, a plane-wave kinetic energy cutoff of 10 hartree and TrouUier-Martins^^ norm- 
conserving pseudopotentials. These calculations have been performed at the theoretical 
lattice constant. To perform the finite difference calculations of the Raman polarizabilitics, 
changes of the electronic dielectric tensor were computed for atoms displaced by ±1% of the 
unit cell parameter along the cartesian directions. 

In case of rhombohedral BaTiOa, we used a 10 x 10 x 10 grid of special k-points, a plane- 
wave kinetic energy cutoff of 45 hartree and extended norm-conserving pseudopotentials^^. 
Since the ferroelectric instability is quite sensitive to the volume of the unit cell and tends 
to disappear due to the volume underestimation of the LDA^^, we chose to work at the ex- 
perimental lattice constants. At the opposite to the lattice parameters, the atomic positions 
have been relaxed : the residual forces on the atoms were smaller than 5 • 10~^ hartree/bohr. 

It was shown by Gonze, Ghosez and Godby^^ that an accurate functional for the exchange- 
correlation energy in extended systems should depend on both the density and the polar- 
ization. The LDA used here neglects this polarization dependence and may consequently 
indroduce significant relative errors when studying the response of a solid to an electric 
field. In case of the second-order derivatives, the LDA usually yields an overestimate of the 
dielectric tensor (as large as 20 % in BaTi03)^^. At the opposite, no clear trends have been 
reported yet concerning non-linear optical properties such as Xifi-^'"^^ 

In LDA calculations, it is common practice to apply a scissors correction^^ to compensate 
the lack of polarization dependence of the exchange-correlation functional. In case of non- 
linear optical properties, such a correction can be applied at different levels. On the one 
hand, we can compute the non-linear optical susceptibilities (Eq. (22)) using a scissors 
operator for the first-order wavefunctions^^ . On the other hand, in the computation of the 
EO coefficients, we can use a scissors corrected refractive index in Eqs (39) and (40). The 
influence of these corrections will be discussed below. 

B. Non-linear optical susceptibilities and Raman polarizabilities of semiconduc- 
tors 

In order to illustrate the computation of third-order energy derivatives, we performed a 
series of calculations on various cubic semiconductors. In these compounds, the non-linear 
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optical susceptibility tensor dijk and the Raman susceptibility tensor aij only have one 
independent element di23 and ai2. Also, instead of ai2 it is customary to report the Raman 
polarizabihty^'^ defined as 

a = a/ iJ,flai2 (46) 

where /i is the reduced mass of the two atoms in the unit cell. 

The formalism of Sec. II involves an integration over the BZ and a differentiation with 
respect to k. In practical calculations, these operations must be performed on a discrete mesh 
of special k-points. As we explained in Sec. II, the discretization can either be performed 
before (PEAD) or after (DAPE) the perturbation expansion of the energy functional Eq. 
(6). Up to know, the applications of the present formalism to real materials^'^° made use of 
the DAPE formula of the third-order energy. The only application of the PEAD formula has 
been reported by Nunes and Gonze^^ on a one-dimensional model system. These authors 
observed that the PEAD formula converges better with respect to the k-point samphng 
than the DAPE formula. In order to compare the performance of these two approaches on a 
realistic case, wc applied both of them to compute the non-linear optical susceptibility di23 
of AlAs. We performed a series of calculations on a n x n x n grid of special k-points. As 
can be seen on Figure 1 the PEAD formula converges much faster than the DAPE formula. 
Therefore, the PEAD formulation has been applied to obtain the results presented below. 
It is the one that is actually available in the ABINIT code. 

40 
35 

> 

E 

B 30 

CO 
CVI 

t5" 

25 
20 

2 4 6 8 10 12 14 16 18 20 
n 

FIG. 1: Non-linear optical susceptibility di23 (pm/V) of AlAs for various grids of n x n x n special 
k-points. 

In Table I, we report the non-linear optical susceptibilities of the cubic semiconductors 
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AlAs and AlP. Our results are in close agreement with the values obtained by Dal Corso 
and co-workers^ who applied the 2n+ 1 theorem within the DAPE formahsm, the results 
of Levine and co-workers^^ who used a "sum over excited states" technique and the values 
obtained by Souza and co-workers^^ who followed a finite electric field approach. The values 
in the lower part of Table I have been obtained using a scissors correction. Our methodology 
provides a correction similar to what is reported by Levine and Allan^°. 

The scissors correction decreases the value of the non-hnear optical susceptibilities in 
agreement with the discussion of Ref. 58. To the authors knowledge, no experimental data 
are available for AlAs and AlP. For other cubic semiconductors, it is however not clear that 
the use of a scissors correction improves the agreement with the experiments^ and will even 
have a negative effect when the LDA underestimates the experimental value. In addition, it 
is not straightforward to isolate the error of the LDA on the non-linear response functions 
from other sources of errors. Other factors have a similar strong influend on xifi the 
scissors correction. For example, the values of the non-linear optical susceptibilities strongly 
depend on the pseudopotential^ or on the error on the unit cell volume^^'^^ that is usually 
underestimated in LDA calculations. 

TABLE I: Non-linear optical susceptibilities di23 (pm/V) of some semiconductors. The values in 
the lower part of the table have been obtained using a scissors (SCI) correction. 



Method AlAs AlP 



2n+ 1 theorem (present) 


35 


21 


2n+ 1 theorem ° 


32 


19 


Finite fields 


32 


19 


Sum over states 


34 


21 


2n+ 1 theorem -|- SCI (present) 


21 


13 


Sum over states -|- SCI 


21 


13 



"Ref. 9 
*Rcf. 57 
"^Ref. 30 



We also computed the Raman polarizabilities of the transverse (TO) and longitudinal 
optical (LO) phonons of various semiconductors. In addition, we performed finite difference 
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calculations of the dielectric tensor with respect to atomic displacements. Our results are 
summarized in Table II where we also report the results of Deinzer and Strauch^° (DS), 
Baroni and Resta^^ (BR) as well as the experimental result of Wagner and Cardona^^ for 
Si. The agreement between our results and those obtained in previous works is quite good. 
In addition, the results we obtained from the 2n + 1 theorem closely agree with the finite 
difference calculations giving us some indication of the numerical accuracy of the implemen- 
tation. 

The Raman polarizabilities of the TO and LO modes are different. As it is discussed 
in Sec. IIIB, this difference is attributed to the macroscopic electric field associated to a 
longitudinal polar lattice vibration. On the one hand, this field modifies the dynamical 
matrix at q ^ 0. The eventual related modification of the eigenvectors of the LO modes 
may imply a first change of the Raman susceptibility. On the other hand, the macroscopic 
electric field itself may induce an additional change of a related to the non-linear optical 
coefficients xiji- the cubic semiconductors, the eigenvectors of the TO and LO modes 
are identical. The difference between the polarizabilities of the TO and LO modes comes 
therefore exclusively from the second term of Eq. (30). 

C. EO tensor in ferroelectric oxides 

In the rombohedral phase of BaTiOa, the EO tensor has four independent elements: ri3, 
r33, r22 and rsi. At the opposite to the dielectric tensor, the EO coefficients can either be 
positive or negative. The sign of these coefficients is often difficult to measure experimentally. 
Moreover, it depends on the choice of the cartesian axes. Experimentally, these axes are 
chosen according to the Standards on Piezoelectric crystals. The z-djds is along the direction 
of the spontaneous polarization and the |/-axis lies in a mirror plane. The z and y axis are 
both piezoelectric. Their positive ends are chosen in the direction that becomes negative 
under compression. The orientation of these axes can easily be found from pure geometrical 
arguments. Unfortunately, these arguments do not allow to determine the direction of the 
?/-axis. Therefore, we apphed the methodology of Ref. [60] to compute the piezoelectric 
tensor from finite differences of the Berry phase polarization. Our results are reported in 
the frame where the piezoelectric coefficient e22 and ess are positive. 

These coefficients, as well as their decomposition on the individual phonon modes and 
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TABLE II: Raman polarizabilities of the transverse (TO) and longitudinal (LO) optical modes 
(A2) of some semiconductors. 

Si AlAs (TO) AlAs (LO) AlP (TO) AlP (LO) 

2n + 1 Theorem 

Present 20.02 
DS 23.56 
Finite differences 

Present 20.17 
DS " 20.44 
BR 26.16 
Experiment 23 =b 4 

''Ref. 10 
''Ref. 37 
^Ref. 59 
''Ref. 63 
'^Ref. 64 

their electronic part, are reported in table IIL All EO coefficients are positive. As it is 
the case for the tetragonal phase^^, the modes that have the strongest overlap with the soft 
mode of the paraelectric phase dominate the amphtude to the EO coefficients. Moreover, 
the electronic contribution is found to be quite small. 

As we discussed in the previous sections, linear and non-linear optical susceptibilities are 
sometimes relatively inaccurate within the LDA. In this context, it is interesting to investi- 
gate the error due to the use of the LDA optical dielectric constants in the transformation 
Eq. (37). Unfortunately, we could not find any experimantal data on the EO coefficients 
in the rhombohedral phase of BaTiOs. In Ref. 46, we studied the EO coefficients of ferro- 
electric LiNbOs and tetragonal BaTiOs and PbTiOa and found an overall good agreement 
between theory and experiment. In Table IV, we report the EO coefficients of these com- 
pounds as well as the values obtained using a scissors corrected optical dielectric constant. 
No scissors correction has been applied for the non-linear optical susceptibilities of these 
compounds that are required to compute the electronic contributions. 

The effect of this correction is more important for the perovskite compounds than for 



8.48 
7.39 



12.48 



4.30 
5.13 



7.46 



8.59 
5.64 



4.25 
4.44 
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TABLE III: Decomposition of the clamped EO tensor (pm/V) in the rhombohedral phase of 
BaTiOa- Reported are the contributions of individual zone-center phonon modes and the elec- 
tronic contribution. The phonon frequencies are reported in cm~^. 





Ai -modes 






E-modes 




CO 


^13 


^33 


CO 


'-22 


r'' 


168 


0.65 


2.16 


163 


0.79 


5.15 


253 


13.82 


27.32 


202 


5.40 


19.16 


509 


1.31 


2.05 


293 
469 


0.01 
0.24 


-0.02 
0.65 


Elect. 


1.15 


2.95 




0.12 


1.24 


Tot. 


16.93 


34.48 




6.56 


26.18 



LiNb03, for which the LDA bandgap and optical dielectric constants are in reasonable agree- 
ment with the experiments^ . For BaTiOa, we tested the optical dielectric tensor obtained 
from the scissors correction that modifies the LDA bandgap to its experimental value^-*^: we 
obtain 1^3 = 12.68 pm/V and rgg = 30.84 pm/V in closer agreement with experimental data. 
However, such an improvement is not a general rule. In PbTiOa, a scissors shift that correct 
the LDA bandgap fails to correct the LDA optical dielectric constant (we obtain en — 5.81 
and £33 = 5.51 while the experimental values are 6.63 and 6.64^^) and yields r^3 = 14.24 
pm/V and r^3 = 8.94 pm/V. Using the experimental dielectric constants, we obtain r^3 = 
10.92 pm/V and r33 = 6.16 pm/V in better agreement with the experiment. 

V. CONCLUSIONS AND PERSPECTIVES 

In this paper, we presented the general framework for the computation of third-order 
energy derivatives within DFT. Our formalism makes use of the 2n + 1 theorem and the 
modern theory of polarization. Focusing on derivatives that are characterized by a zero 
wavevector and that involve either three electric fields or two electric fields and one atomic 

displacement, we described the computation of non-linear optical susceptibilities, of Raman 
scattering efficiencies of TO and LO phonons and of the EO tensor. 

The computation of the Berry phase polarization involves a derivative of the wavefunc- 
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TABLE IV: Effect of a scissors correction on tlie EO coefficients of LiNbOs (clamped and un- 
damped cases), and the tetragonal phases of PbTiOs and BaTiOs (clamped cases only). The 
dielectric tensor required to perform the transformation Eq. (37) has been computed whitin the 
LDA {elda) and using the LDA with a scissors correction (ssci)- No scissors correction has been 
used to compute the non-linear optical susceptibilities that determine the electronic contribution 



to the EO coefficients. In case of PbTiOs, we also use the experimental dielectric tensor {£exp) to 
compute the EO coefficients. The values are compared to the experimental results. 






V \ 'J 

^ 1.) 




J'99 


1'- 1 


T.iNbOo 


^LDA 


Q fi7 








(clamped) 


£SCI 


10.37 


28.89 


4.88 


16.02 




iliXp. 


O.D 


on o 
oO.o 






LiNbOg 


£lda 


10.47 


27.08 


7.53 


28.61 


(undamped) 


Ssci 


11.23 


29.06 


8.08 


30.69 




Exp.72 


10 


32.2 


6.8 


32.6 




Exp.^3 






9.89 




PbTiOa 


£LDA 


8.98 


5.88 




30.53 


(clamped) 




14.24 


8.94 




47.39 




Sexp 


10.92 


6.16 




34.45 




Exp.74 


13.8 


5.9 






BaTiOs 
(damped) 


£LDA 
£SCI 


8.91 
12.68 


22.27 
30.84 







Exp.^^ 10.2 40.6 

Exp.^^ 8 28 



tions with respect to their wavevector. In practice, this differentiation is computed on a grid 
of special k-points. The perturbation expansion can either be performed before (DAPE) or 
after (PEAD) the discretization, leading to two mathematically distinct expression of the 
third-order energies. We used both of them to compute the non-linear optical susceptibility 
of AlAs and we have shown that the PEAD formulation converges faster with respect to the 
k-point sampling. 
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We have computed the non-hnear optical susceptibihties and Raman polarizabihties of 
some cubic semiconductors as well as the EO tensor in the rhombohedral phase of BaTiOs. 

Finally, we have studied the effect of a scissors correction on the EO coefficients and the 
non-linear optical susceptibilities. At the opposite to the dielectric tensor, we did not find 
a systematic improvement of the results by using this correction. 

We can figure out several applications of the methodology presented in this work. Com- 
bined with the calculation of phonon frequencies and infrared intensities, the computation 
of Raman efficiencies can be a useful complementary tool for the interpretation of experi- 
mental spectra. Furthermore, the computation of the EO tensor from first-principles can 
guide the tuning of the EO properties and help designing new efficient EO materials. This 
could reveal particularly helpful since accurate optical measurements require high quality 
single crystals not always directly accessible. 
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APPENDIX A: EXPRESSIONS OF THE CLAMPED AND UNCLAMPED EO 
TENSORS 

1. Macroscopic approach 

As it is discussed in Sec. HID, the optical properties of a compound are modified by an 
electric field or a mechanical constraint (a stress a^j^i, or a homogeneous strain 77^,/). At 
the linear order, the variations of ej^ can be described using either the variables (£-y, 'q^^) or 
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3 3 

7=1 /i,i'=l 
3 3 

7=1 /i,!^=l 

where r^^vy and rfj^ are resp. the clamped (strain free) and undamped (stress free) EO 
coefficients, Pij^u are the elasto-optic (strain-optic) coefficients and Trj^^i/ are the piezo-optical 
(stress-optical) coefficients. In order to relate Eqs. (Ala) and (Alb), we can express the 
strain as beeing induced by the stress or by the electric field (converse piezoelectric effect) 

3 3 
/i',i/'=l 7=1 

where -S'^j/^v' are the elastic compliances and d-y^,^ the piezoelectric strain coefficients. 

If we assume for example that the unit cell is free to relax within the electric field (stress- 
free mechanical boundary conditions) we can either use Eq. (Alb) (in which case the second 
term of the right-hand-side is zero) or Eq. (Ala) to compute A{e~^)ij. In the latter case, 
the strain induced by the electric field can be obtained from the second term of the right- 
hand-side of Eq. (A2) 

3 
7=1 

3 3 3 

~ Z^'"ij'7^7+ Z^ 'y^^Pijfj.ud^f^u^-/- (A3) 
7=1 M>'^=1 7=1 

Using this identity, we obtain the following relation between the undamped and the clamped 
EO coefficients 

3 

M,i^=l 

2. Microscopic approach 

In order to derive the expressions of the clamped and undamped EO tensor of Sec. HID, 
we use a Taylor expansion of the electric enthalpy^'' F. Similar developments have already 
been applied to determine the lattice contribution of the static dielectric tensor and of the 
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piezoelectric tensor^^'^^. They are based on an expansion of F up to the second order in the 
atomic coordinates R^a-i the homogeneous strain 77^,^ and the macroscopic electric field 
In this section, we extend these developments to the third order. 

The electric enthalpy of a solid in an electric field is obtained by the minimization 

F(£) = min F (R, r/, 5) . (A5) 

We denote R(^), ^(^) the atomic positions and the strain that minimize F at constant 
£ and Rq, ryo (= 0) their values at £ = 0. For small fields, we can expand the function 
F (R, ?7, £) in powers of £ around £ = 0: 

F (R, r],E) = F (R, 77, 0)-O ^ (R, r,) 8,-^ J] (R, r/) 8,£j~ J] xgi (R, v) Si8^Sk+- 

1=1 jj=l i,j,k=l 

(A6) 

where Q is the volume of the primitive unit cell in real space and P(R, 77), £jj 
X^^l (R, v) the macroscopic polarization, electronic dielectric tensor and non-linear op- 
tical coefficients at zero macroscopic electric field and for a given configuration (R, 77). At 
non-zero field, these quantities arc defined as partial derivatives of F with respect to £. For 
example, the electric field dependent electronic dielectric tensor can be computed from the 
expression 

(A7) 



An d'^F 

eij{R{8),ri{8),8)^ 



88 id 8 j 



R(£),7?(f),f 



Let Ti^a — R-KQ ~ R-o./ta the displacement of atom k along direction a and r^^ (resp. 77^^) 
the first-order modification of the atomic position (resp. strain) induced by a perturbation 
A 



A=o dX 



(AJ 

A=0 



In the discussion that follows, we will study the effect of an electric field perturbation and 
a strain perturbation on the electric enthalpy F in order to obtain the formulas to compute 
the elasto-optic coefficients as well as the clamped and the undamped EO tensors. 
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a. Elasto-optic coefficients (£ = 0) 



The elasto-optic tensor can be computed from the total derivative of the dielectric tensor 
with respect to r^^^ at zero electric field 



dSij (R, T], 0) 



dsij (R, 77) 



IJ,V 



Ro,m 



(A9) 



The derivative in the first term of the right hand side is computed considering the ionic 
cores as artificially clamped at their equilibrum positions. The remaining terms represent 
the ionic contribution to the elasto-optic tensor. They involve derivatives of the linear 
dielectric susceptibility xfj with respect to the atomic positions that have to be multiplied 
by the first-order strain induced atomic displacements rHa [Eq. (A8)]. To compute these 
quantities we use the fact that F is minimum at the equilibrum for an imposed strain 77. 
This condition implies 

^^('^••'' =0. (AlO) 



R(7?),7? 



Since we are interested in first-order atomic displacements we can write Tkq; (77) = 
TKa^r/^i/ -I- 0{rf'). Solving the extremum equation (AlO) to the linear order in 77, 
we obtain 



92F(R,77) 



d'F{R,r)) 



(All) 



The second derivatives on the left side of Eq. (All) define the matrix of interatomic force 
constants at zero macroscopic electric field which enables the computation of the transverse 
phonon frequencies uj^ and eigendisplacements ^^(^0;). By decomposing rHa in the basis 
of the zone-center phonon-mode eigendisplacements 



E 



(A12) 



and using Eqs. (24), (25) we derive the following expression for the first-order strain induced 
atomic displacements 



-1 d^F(R,r]) 



where 



d'^F (R, rj) 



d^F (R, 77) 



= E 



(A13) 



(A14) 
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If we introduce Eqs. (A12) and (A13) into Eq. (A9) and use the definition of the Raman 
susceptibihty Eq. (28) and the transformation Eq. (37), we finally obtain the formula to 
compute the elastooptic tensor 

■1 deij{R,r]) 



■^'^^^ njn] dfj^^ 



47r ^a-92F(R,r^) ^^^^^ 



Roflo 



To simplify, we write Eq. (A15) in the principal axes of the crystal under investigation. A 
more general expression can be obtained from Eq. (37). 

Eq. (A15) is different from the approach used previously by Dctraux and Gonze to study 
the elasto-optic tensor in a-quartz™. The authors of Ref. 70 used finite differences with 
respect to strains to compute the the total derivative of Sij. In their approach, the atoms 
where relaxed to their equilibrum positions in the strained configurations. In case of Eq. 
(A15), the first term of the right-hand-side is computed at clamped atomic positions while 
the effect of the strain-induced atomic relaxations is taken into account by the second term. 



b. Clamped EO coefficients (rj = 0) 

The clamped EO tensor can be computed from the total derivative of the electric field 
dependent dielectric tensor Eq. (A7) with respect to £ 



dSij (R, r]o,S) 



dSij {Ro,r)o,£) 



Ro,£=0 



e=o 



(A16) 



-Ro 



The derivative in the first term is computed considering the ionic cores as artificially clamped 
at their equilibrum positions. This term represents the bare electronic contribution to the 
EO tensor that can be computed from the non-linear optical coefficients 

dsij {'Ro,r]o,£) 



A;=7 



(A17) 



related to a third-order partial derivative of F 



(2) ^ (2) X ^ -l d^F{Ro,rio,£) 



(A18) 



£•=0 



The remaining terms in Eq. (A16) represent the ionic contribution to the EO tensor. They 
involve derivatives of the linear dielectric susceptibility xi]^ with respect to the atomic posi- 



tions that have to be multiplied by the first-order electric field induced atomic displacements 
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I Ka 



[Eq. (A8)]. To obtain these quantities, we proceed the same way as in case of the elasto- 



optic tensor. Using the equihbrum condition 



dr, 



1=1 



Stt ^ 



a£ij (R, rjo) 



9Tk 



and expanding r^a to the first-order in the electric field, we obtain 



E 



92F(R.//o.O) 



Ro 



(A19) 



(A20) 



Ro 



This expression is similar to Eq. (All). The second-order derivatives of F on the left side are 
the interatomic force constatnts and the derivative of the zero field polarization with respect 
to r^a on the right side is the Born effective charge tensor Z*^^^ of atom k. Decomposing TkZ 
in the basis of the zone-center phonon-mode eigendisplacements [Eq. (A12)] and using the 
orthononormality constraint Eq. (25) we derive the following expression for the first-order 
electric field induced atomic displacements 



(A21) 



If we introduce Eqs. (A17) and (A21) into Eq. (A16) we finally obtain the formula to 
compute the total derivative of the dielectric tensor 



deij (R, £) 



Ro,S=0 



1 



A;=7 



9x!j'(R) 



(A22) 

Using the definition of the Raman susceptibility [Eq. (28)], the mode polarity [Eq. (41)] 
and the transformation [Eq. (37)] we obtain the expression of the clamped EO tensor 



U7 



-StT (2) 



An 



1=1 



E 



out, 



(A23) 



As in case of the elasto-optic tensor [Eq. (A15)], we have written Eq. (A23) in the principal 
axes of the crystal under investigation. 



c. Undamped EO tensor (a = 0) 

In order to compute the undamped EO tensor, we have to take into account both the 
electric field induced atomic displacments rfa and the electric field induced strain rj^Z when 
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computing the total derivative of e. 



^3 



dSij (Ro,r7o,£^) 



i?O,??O,f=0 



3 

+47r 



ax!," (Ro,»)) 



(A24) 



The electronic contribution [first term of Eq. (A24)] is the same as for the clamped EO 
tensor. It can be computed from the non-linear optical coefficients [Eq. (A17)]. To compute 
rfo and rifZ-, we can use an equilibrum condition similar to Eq. (A19) where we require that 
the first-order derivatives of F with respect to t^q and 77^1^ vanish. Expanding t^q and 77^^,/ 
to the first-order in the electric field, we obtain the system of coupled equations [see also 
Ref. 71] 



E 



92F(R,r7,0) 



K'a' o 



R-0,% 



^^7 _ Q 



(A25a) 



R-o,»?o 



/if 



(A25b) 



Because of the coupling between rfa and -qfn^ defined by the mixed second-order derivatives 
Q^^^^ , the second term of the right hand side of Eq. (A24) is different from that of Eq. 
(A16). That means that the sum of the first and second term of Eq. (A24) is not identical 
to the clamped EO coefficients r^-^. Moreover, the third term of Eq. (A24) is different from 
the piezoelectric contribution of Sec. A 1. 

In order to obtain the decomposition of r^-^ into electronic, ionic and piezoelectric contri- 
butions defined previously, we can solve Eq. (A25a) for T^a- In the basis of the zone-center 
phonon mode eigendisplacements we can write 

Pn,^ ^_^^^F{Y^,7],Q) 



n 



R.o,»?o 



,£7 



(A26) 



If we insert this relation into Eq. (A24) and use the transformation Eq. (37) we obtain the 
following expression of the undamped EO tensor in the principal axes 



-Stt 



(2) 



An 



"I '"J 



i=7 



E 



E 



dn 



Ro,'no,£=o 



d'^F (R, 77, 0) 



dTmdr]^ 



Ro,r)o,£=0 
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The sum of the first and second term of the right-hand side of Eq. (A27) is equal to the 
clamped EO coefficient r^^^. The product of the conversion factor times the bracket in the 
third term of Eq. (A27) is equal to the elasto-optic coefficient pij^i, [Eq. (A15)]. Finally, by 
definition of the converse piezoelectric effect, rj^l is equal to the piezoelectric strain coefficient 
d'^^lu■ We obtain thus the following expression of the undamped EO coefficients that is equal 
to the one derived in Sec. A 1 from pure macroscopic arguments 



It is worth noticing that the so-called piezoelectric contribution does not only take into 
account the change of the linear optical susceptibility with strain (third term of the right 
hand side of Eq. (A24)) but also includes the modification of the ionic contribution, with 
respect to the clamped case, that is associated to the modification of the ionic relaxation 
induced by the strain. 
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